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We explain how to implement, in the context of projected entangled-pair states (PEPS), the 
general procedure of fermionization of a tensor network introduced in [P. Corboz, G. Vidal, Phys. 
Rev. B 80, 165129 (2009)]. The resulting fermionic PEPS, similar to previous proposals, can be 
used to study the ground state of interacting fermions on a two-dimensional lattice. As in the 
bosonic case, the cost of simulations depends on the amount of entanglement in the ground state 
and not directly on the strength of interactions. The present formulation of fermionic PEPS leads to 
a straightforward numerical implementation that allowed us to recycle much of the code for bosonic 
PEPS. We demonstrate that fermionic PEPS are a useful variational ansatz for interacting fermion 
systems by computing approximations to the ground state of several models on an infinite lattice. 
For a model of interacting spinless fermions, ground state energies lower than Hartree-Fock results 
are obtained, shifting the boundary between the metal and charge-density wave phases. For the 
t — J model, energies comparable with those of a specialized Gutzwiller-projected ansatz are also 
obtained. 
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I. INTRODUCTION 

Strongly correlated fermionic systems, responsible 
for relevant many-body phenomena such as high- 
temperature superconductivity, the fractional quantum 
Hall effect or metal-insulator transitions, represent one of 
the most important theoretical challenges in condensed 
matter physics. Among the simplest possible models 
of interacting fermions in a 2D lattice is the Hubbard 
model,'ij which is believed to be one of the keys to un- 
derstanding the theoretical riddle of high-temperature 
superconductivity,'^ and which serves as a good exam- 
ple to illustrate the nature and scale of the difficulties 
encountered. In spite of a titanic effort by the condensed 
matter community spanning several decades, still today 
the phase diagram of the 2D Hubbard model and its rela- 
tion to high-temperature superconductors remain highly 
controversial. 

In the absence of exactly solvable models, accurate 
numerical simulations are essential in order to gain fur- 
ther insight into the physics of strongly correlated sys- 
tems. While quantum Monte Carlo (QMC) techniques 
are very powerful in simulating bosonic systems, they suf- 
fer from the so-called negative signvroblem in the case 
of fermionic and frustrated models.f^ On the other hand, 
generic ID lattice systems can be accurately addressed 
with the density matrix renormalization group (DMRG) 
method,!^ but this approach scales inefficiently with the 
lattice size in 2D systems. Recent progress in the simula- 
tion of 2D fermionic models has been made with a variety 
of methods.'SlS] However, results obtained with different 
methods are often inconsistent, highlighting the need for 
further improvement and for alternative approaches. 

A promising new route to studying strongly correlated 
fermion systems in a 2D lattice, presently under intense 



investigation,'2Hill is based on using a tensor network as 
ground state variational ansatz. For bosonic (e.g. spin) 
2D lattice models, tensor network ansatze include pro- 
jected entangled-pair states (PEPS) for inhomogeneou^^ 
and homogeneous systems ,^^"221 and the multi-scale en- 
tanglement renormalization ansatJ^ (MERA). [Homoge- 
neous PEPS are also known with names such as (vertex) 
tensor product state^^^^\. The interest in these ap- 
proaches resides in the fact that they manage to retain 
some of the useful features of DMRG and QMC, while 
avoiding their main shortcomings. Indeed, PEPS and 
MERA approaches are free from the negative sign prob- 
lem that prevents the application of QMC to fermionic 
and frustrated models. At the same time, and unlike 
DMRG, both PEPS and MERA can efficiently represent 
ground states of 2D lattice models. In addition, com- 
pared to other variational approaches, PEPS and MERA 
are relatively unbiased towards specific ground states. 
Still at an early stage of development, the major lim- 
itation of these methods is that the cost of simulations 
increases sharply with the amount of entanglement in the 
ground state. This limits the range of models that can 
be analyzed accurately at present. Nevertheless, several 
systems of frustrated antiferromagnets beyond the reach 
of DMRG and QMC have already been addressed-l^lHlIl 

In recent months, generalizations of tensor network al- 
gorithms to fermionic systems have been put forward in- 
dependently by several groups .I^Hn] As a result, it is now 
possible to study interacting fer mions i n 2D lattic es both 
within the context of the MER^aEMl and PEPS.ESMIHI 
The fundamental new step, common in all the proposals, 
is to incorporate the fermionic character of the ground 
state wave function directly into the ansatz. This is ac- 
complished by considering a network of fermionic opera- 
tors, that is, a set of linear maps, made of anticommuting 
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operators, that are connected according to a network pat- 
tern, as first proposed in Refs. i9|lli for the MERA and 
in Ref.[inifor the PEPS. 

In actual computations, one is still forced to store and 
manipulate tensors (i.e. multi-dimensional arrays of co- 
efficients) corresponding e.g. to the expansion coeffi- 
cients of the fermionic maps. The process of ^fermion- 
ization' of a tensor network algorithm, i.e. its extension 
to fermionic systems, can in practice be achieved and 
visualized in a variety of (ultimately equivalent) ways, 
depending on how the underlying network of fermionic 
operators is translated into a set of tensors and rules 
for their manipulation. Exa mples include the use of a 
Jordan- Wigner transformation,I^Ei]or the introduction of 
additional bond indices between tensors.^iSl A particularly 
simple form of ^ fermionization^ of tensor networks was in- 
troduced in Ref . [HI where it was applied in the context of 
the MERA. We emphasize that Ref. is based on refor- 
mulating previous work by Corboz, Evenbly, Verstraete 
and Vidal on fermionic MERA,!^ which in turn had its 
origins in a key observation by Verstraete.'^Sl 

The fermionization procedure of Ref. 1121 which applies 
to any tensor network, is remarkedly simple. It does 
not require the introduction of a Jordan- Wigner trans- 
formation in the bond indices, or to have to explicitly 
keep track and dynamically modify a global fermionic 
order; neither does it require the introduction of addi- 
tional bond indices in the tensors. Instead, the fermionic 
character of the tensor network is engraved in two sim- 
ple rules: (i) use of parity invariant tensors and (ii) re- 
placement of line crossings with so-called fermionic swap 
gates. The net result of applying these rules is a modified 
variational ansatz that can be manipulated using stan- 
dard tensor network operations (tensor multiplications, 
etc.), thus producing a straightforward fermionic version 
of existing tensor network algorithms. Importantly, the 
computational cost of bosonic and fermionic algorithms 
scales in the same way with the amount of entanglement 
in the ground state.'^ This remarkable result was also 
independently derived in Ref. [T31 

This paper has two main goals. The first is to explain 
how to obtain PEPS algorithms for fermionic systems 
by applying the above ^fermionization^ rules to existing 
bosonic PEPS algorithms. Fermionic PEPS were origi- 
nally proposed by Kraus, Schuch, Verstraete and Cirac in 
Ref. Wand have also been discussed by Barthel, Pineda 
and Eisert in Ref. 1131 Our formulation of fermionic 
PEPS must be, at some level, equivalent to those pro- 
posals. However, the present formulation, which is based 
on previous independent work,I^EIIis remarkably straight- 
forward and appears to be comparatively much simpler. 
In particular, it allowed us to numerically implement a 
fermionic PEPS algorithm for infinite systems by only 
introducing a small number of changes to existing code 
for bosonic systems. 

A second main goal of this paper is to demonstrate 
the usefulness of fermionic PEPS. In spite of the several 
existing formulations of fermionic PEPS,^" ^and with 



the exception of Ref. [TU where some qualitative results 
are reported for the t — J model, no evidence has been 
presented yet showing that fermionic PEPS are a good 
variational ansatz for interacting fermion systems. [No- 
tice, however, that Ref.[TU]shows that Gaussian fermionic 
PEPS can represent states of non-interacting fermions]. 
Here we do present such evidence, in the form of ground 
state computations for several 2D models. 

Specifically, we use a f ermio nic version of the infinite 
PEPS (iPEPS) algorithnil^iao] ^ address models on an 
infinite square lattice. First, results for free spinless 
fermions are compared with the corresponding exact so- 
lution, showing that a PEPS with small bond dimension 
is capable of reproducing the ground state energy with 
several digits of accuracy. Then a model of interacting 
spinless fermions is addressed. Qualitatively, the sim- 
ulation reproduces the phase diagram predicted within 
Hartree-Fock, with metal and charge-density wave phases 
separated by a line of first order phase transitions. At a 
quantitative level, however, we obtain ground state en- 
ergies that are lower than those obtained with Hartree- 
Fock, and this shifts the boundary between phases sig- 
nificantly. Finally, for the t — J model, we obtain ground 
state energies that are close to those of a specialized 
Gutziller-projected ansatz. 

The rest of the paper is organized as follows: Sec. 
reviews the PEPS formalism for bosonic systems and the 
general fermionization procedure of tensor networks in- 
troduced in Ref. 112] which is then applied to PEPS al- 
gorithms. Sec. |III| considers in more detail the fermionic 
version of the iPEPS algorithm for infinite 2D lattices, 
which was employed to obtain the benchmark results pre- 
sented in this paper. Sec. |IV| describes ground state cal- 
culations for systems of free and interacting fermions in 
an infinite 2D lattice. Sec. [V] contains some conclusions, 
while Appendix [A] defines generalized fermionic opera- 
tors and Appendix [B] describes in detail one step of the 
update in the fermionic iPEPS algorithm. 

Note on terminology. — For the purposes of this paper, 
a tensor is simply a multi-dimensional array of complex 
coefficients, and a tensor network is a set of tensors some 
of whose indices are connected according to a network 
pattern, where being connected means that there is a 
sum or trace over that index, in the sense of tensor multi- 
plication. Accordingly, in this paper a bosonic/fermionic 
tensor network is a tensor network used in the context 
of simulating a bosonic/fermionic system. Thus, in the 
present formulation a fermionic PEPS is simply a "tensor 
network that serves as a variational ansatz for fermionic 
systems" . It is different from a bosonic PEPS in the pres- 
ence of special gates called fermionic swap gates (and in 
that its tensors are necessarily parity preserving). In par- 
ticular, even though the rules used to create a fermionic 
tensor network, as introduced in Ref. [T^ and reviewed 
here, were obtained by studying how to mimic a network 
of fermionic operators (that is, of operators that obey 
anticommuting relations), here a fermionic PEPS is not 
a network of fermionic operators. One of the merits of 



the present formulation is precisely that it replaces the 
considerable complexity involved in dealing with a net- 
work of fermionic operators with a simple set of rules. 
In particular, it avoids having to explicitly define, keep 
track and dynamically modify a fermionic order for the 
bond indices. The equivalence between our formulation 
of a fermionic tensor network and a network of fermionic 
operators was already established in Ref . il2, for the case 
of the MERA. A general derivation of this equivalence 
would distract from the purpose of this paper and will 
be presented elsewhere. 

II. FERMIONIZATION OF PEPS 

The goal of this section is to introduce a fermionic ver- 
sion of bosonic PEPS algorithms,l2^^so that they can be 
applied to simulate fermionic systems in a 2D lattice. We 
start by reviewing some key aspects of the PEPS formal- 
ism for bosonic systems. This allows us to introduce the 
notation and the diagrammatic representation of tensors 
used throughout Sees. [ll] and |III[ Then we describe the 
fermionization rules of Ref. 1121 which we also extensively 
review. We apply these rules to obtain a fermionic PEPS 
ansatz (see also Refs. ll0ll3ll4p . and provide a discussion 
of how fermionic PEPS algorithms can be obtained by 
modifying existing bosonic PEPS algorithms. 

A. Bosonic lattice system 

Let us consider a quantum many-body system in a 
lattice £ made of N sites, labelled by an integer k S 
{1, 2, • • • , N}. Each site k £ £ is described by a complex 
vector space V of finite dimension d, with basis states 
{|s)}s=i^... ^d- The vector space V could represent, for in- 
stance, the possible states of a quantum spin sitting on 
that site of £. The system is further characterized by a 
local {bosonic) Hamiltonian. This is a Hermitian opera- 
tor H : V*^^ — > Y'g'^ i^jiat (when expressed in terms of 
bosonic operators, i.e. operators that commute when act- 
ing on different sites) decomposes as a sum of terms each 
involving only a small number of sites. Let \'^) € 
be a pure state, 

l*>= ^si.....s„|siS2---Sn), (1) 

where index Sk labels a basis on site k £ C. 

A task of interest is to compute a specific state j^*) 
somehow related to H, e.g. its ground state, and to 
evaluate the expectation value (^'|6|5') of some local ob- 
servable o. However, representing a vector \'i>) £ Y'^'^ 
requires a number of complex coefficients '^siS2---sisi that 
grows exponentially in N. This poses a serious com- 
putational challenge. Exact diagonalization techniques 
are only affordable for small systems (e.g., at most N « 
30 — 40 for d = 2), and alternative numerical strategies 
are required to analyze large systems. 
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FIG. 1: (Color online) (a) Diagrammatic representation of the 
tensor 'I' with coefficients ^'^ia' -is for a state l^*) G V®^ of a 
3x3 square lattice jC. This tensor is expressed in terms of a 
PEPS made of a set of 9 tensors {A^^}, one for each site f £ C. 
(b) Bulk tensor A^^ with components j4[^^^^. Notice that the 
legs corresponding to indices u, I, s, d and r emerge from 
the circle in anti-clockwise order, (c) Hermitian conjugate 
A^^'' of the PEPS tensor A^^, represented as its mirror image. 
Notice that the legs corresponding to indices r, d, s, I and 
u of (A^^^)rd3iu emerge again in anti-clockwise order, (d) 
Hermitian conjugates of and its PEPS representation. 



B. Projected Entangled-Pair States 

Projected entangled-pair stateJi^Ell (PEPS) were in- 
troduced as a means to obtain an efficient description 
for some states \^) € V*^ of a 2D lattice £. For con- 
creteness, in this work we consider the case of a square 
lattice, although all the discussions can be extended to 
other type of lattices. To each site k £ £ there corre- 
sponds a vector of integers r = {x{k),y{k)), and we also 
write f£ £ to denote a site of lattice £. 

A PEPS is made of a collection of iV tensors {^I*^}, one 
for each site r £ £, connected through bond indices that 
follow the pattern of links of the lattice £. Upon tracing 
over all bond indices, a PEPS yields a tensor 5* with the 
d^ complex coefficients *siS2 -s« of a state l^*) £ V*^. 

Throughout this paper a diagrammatic representation 
of tensors and tensor networks is used, see Fig. [T] Each 
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tensor is depicted as a shape (circle, square, diamond, 
etc) and its indices as emerging lines. A line connecting 
two shapes (or starting and ending at the same shape) 
denotes an index over which a trace is taken. As an exam- 
ple, Fig. [ija) represents, in the case of a lattice £ made 
of 3 X 3 sites, a tensor with 9 indices, corresponding 
to the coefficients ^'s^s^.^jg of I*) e , followed by a 
PEPS made of 9 tensors {A^^}. 

The number of indices in a tensor A^^ depends on the 
number of nearest neighbors of site f € C, with tensors in 
the bulk having more indices than tensors at a boundary. 



Specifically, a bulk tensor has components A 



■_r\ 
ulsdr ' 



with 



one physical index s and four bond indices u,l,d, r. The 
physical index s labels the basis of the vector space V for 
site f & C and therefore takes d different values, whereas 
each bond index connects the tensor with a tensor in a 
nearest neighbor site and ranges from 1 to D, where D is 
the so-called bond dimension of the PEPS. Correspond- 
ingly, a bulk PEPS tensor is represented by a circle with 
five legs, see Fig. [l|b). As in the rest of the paper, here 
we follow the prescription that the indices of a tensor are 
drawn in anti-clockwise order. Notice also that the open 
indices of the PEPS in Fig. [ija) reach the exterior of the 
tensor network in exactly the same (anti-clockwise) or- 



der that they appear in '^i-^^i 



These notational and 



diagrammatical details are somewhat superfluous in the 
bosonic case (since one can change the order of indices 
in a tensor by simply permuting its components) but will 
become important in the extension to fermions. 

A PEPS on a L X L lattice, where N — L'^, con- 
tains 0{N) bulk tensors, each depending on dD'^ com- 
plex coefficients. Therefore the PEPS is characterized by 
0{NdD^) parameters. If D has a fixed value indepen- 
dent of iV, then the PEPS is indeed an efficient encoding 
of some states \^) G V*^^, and it can be used e.g. as 
a variational ansatz for the ground state of H . How- 
ever, for the PEPS to be a useful ansatz, we also need to 
provide an efficient strategy to optimize its tensors and 
manipulate them in order to extract physically relevant 
information, as discussed next. 



C. Optimization and expectation values 

Two major tasks to be accomplished with a PEPS 
algorithnJISlEni q^-q- ^j) optimization of its 0{NdD^) pa- 
rameters so as to obtain a good approximation to e.g. 
the ground state of the local Hamiltonian H] (ii) given a 
PEPS for a state \^) , computation of expectation values 
(^'|6|'I') of local observables 6, as given by 



mom 



(2) 



These two tasks happen to be closely related. They in- 
volve taking the trace over all the bond indices of a com- 
posite tensor network, an operation referred to as con- 
tracting the tensor network. 
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FIG. 2: (Color online) (a) Scalar product ('I'l*) written in 
terms of tensors and 'I'^. (b) The same scalar product, 
but written in terms of a PEPS together with its Hermitian 
conjugate, (c)-(d) Using the jump move of Fig.lsl this tensor 
network can be modified so that each tensor A^^ is drawn 
next to its Hermitian conjugate A'*^. (e) Tensor network £ in 
terms of reduced tensors a'*^. (f) a'*^ defined in terms of A^^ 
and ^["^t^ Eq. ^. 



An emblematic example of tensor network contraction 
required in a PEPS algorithm concerns the computation 
of the scalar product We need to introduce a few 

more definitions and notation. Let denote the Hermi- 
tian conjugate of a tensor T, obtained by reversing the 
order of the indices of T and taking the complex con- 
jugate of each of its coefficients (and diagrammatically 
represented as the mirror image of T). For instance, ten- 
sor "if^ in Fig. [ijd), corresponding to (^|, 



has coefficients 



W2 



■in\ 



mh 



(3) 



(4) 



Fig. [ijc) represents the Hermitian conjugate A^t of a 
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FIG. 3: (Color online) Jump move: The graphical representa- 
tion of a tensor network is not unique. In particular, a line can 
be dragged over a circle without changing the tensor network 
that is represented. This property, trivial in tensor networks 
for bosonic systems, will have a less obvious analogue in the 
fermionic case. 



bulk tensor A^^, which has coefRcients 

{A^'^'Usiu^A^l,^ (5) 

We build the Hermitian conjugate of the PEPS for ^ 
as the set of tensors {A^^'''} connected according to the 
mirror image of the network of connections in the PEPS 
for see Fig. [ijd). For each site r € £, let us define 
the reduced tensor a^^ in terms of tensor A^"^ and A'^^^ 
by tracing over the physical index s. For instance, in the 
bulk, the reduced tensor a^^ has components 

d 

°'uullddfr = X! ^ulsdriA^'^^)fdslu^ (6) 
s = l 

see Fig.[2jf). Then the scalar product (^|^) results from 
contracting a tensor network £ that consists of all re- 
duced tensors a^^ connected according to the links in C, 
see Fig. [2|^a)-(e), where the (trivial) jump move of Fig. [s] 
is used. 

Contracting the tensor network £ to obtain the scalar 
product (^'1^') comes with a cost that grows exponen- 
tially in the linear size L of lattice C, and therefore cannot 
be accomplished efficiently. A key ingredient of PEPS al- 
gorithms is precisely a strategy to efficiently but approx- 
imately contract the tensor network £, thus producing 
an approximation to (\['|^'). This can be done in several 
ways, depending on the size and topology of lattice C. In 
a finite lattice with open boundary conditions, one can 
use matrix product state (MPS) techniques.!^ In the case 
of a torus, coarse-graining techniques know n as tensor 
entanglement renormalization group (TERG)^'^can be 
used. Finally, in an infinite lattice , both infinite MP^^^ 
and corner transfer matrix (CTM)2212ZI techniques have 
been employed. 

In order to optimize a PEPS so as to approximate the 
ground state of H, as well as to evaluate the expectation 
value of local observables, it is useful to contract cer- 
tain class of tensor networks called environments. The 
environment £^^ of a site r & C is the tensor network ob- 
tained from £ by removing tensor a^'^, and can be used 
to compute the expectation value of a local observable 
acting on that site. 




FIG. 4: (Color online) (a) Tensor network £, made of reduced 
tensors a^^, corresponding to the scalar product in an 

L X L lattice £. with open boundary conditions, (b) Envi- 
ronment £;[''i''2l for two contiguous sites fi,f2 £ £, obtained 
from £ by removing tensors a''^^' and a'""^'. (c) Approximate 
environment consisting of six tensors {Ea}. (d) Tensor 

o of coefficients for a two-site local operator 6, Eq. Q. (e) 
Tensor network made of the approximate environment c;''^^'^' 
and tensors A^^^\ A^^^^\ A^^^\ A^'^^lt .^^^ ^ j^s contraction 
produces an approximation to (^|6|>I'). 



Similarly, a two-site environment £^^^^^^ is the tensor 
network obtained by removing tensors a^^^^ and a^'"^! from 
f , and can be used e.g. to compute the expectation value 
of a two-site observable 

0= Oi2ilJU2b'lJ2>(il»2|, (7) 

acting on sites ri,r2 G C. Figure |4]^b) shows the en- 
vironment £^[''i''2l corresponding to two nearest neigh- 
bor sites ri,r2 S ^- Again, the exact contraction of 
the environment cannot be performed efficiently, but ef- 
ficient schemes, analogous to those employed to con- 
tract £, can be used in order to approximately contract 
£lri,r2]^ The whole environment is in this way approxi- 
mated by a smaller tensor network Q^^^"^^^ made of 6 ten- 
sors {El, ■■■ , Eq}, see Fig. [4jc). These tensors can then 
be connected to tensors yl™, A^^^\ A^^lt and o 

to produce an approximation to (^lol^/), see Fig. |4]je). 

In a system with a Hamiltonian H made of two-site in- 
teractions between nearest neighbors, the approx- 
imate environment Q^^^^^^ is particularly useful. On the 
one hand, it is employed in the optimization of tensors 
and A^^^^ after a gate has been applied on the two 
sites, as part of simulating an imaginary-time evolution 
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according to -ff , 



~ lim 



o-tH 



*0> 



(8) 



which offers one way of obtaining a PEPS approxima- 
tion to the ground state I'i'cs) of H. On the other hand, 
glr-ir2] pg^Q oigQ |-,g yggjj compute the expectation value 

of the energy on that Hnk, {'$\h'^^^^^^ |^), as part of an al- 
gorithm to optimize the PEPS by minimizing its energy, 



mm — 



(9) 



which is another way of obtaining a PEPS approximation 
to the ground state |*gs) of H. We refer to Refs. ll5)17)20l 
for more details. Other PEPS algorithms'^ '^bypass the 
computation of environments. 

This concludes our short review of PEPS algorithms 
for bosonic 2D lattice models. 



D. Fermionic lattice systems 

Let us now consider a fermionic lattice system. For the 
sake of simplicity, we first assume that each site fc € £ is 
described by a complex vector space V of dimension d = 2 
that is associated to a fermionic annihilation operator , 
with anticommutation relations 



cjjCfc' + Ck'cl = 6kk' 
CkCk' + Ck'Ck = 0, 



(10) 
(11) 



[we will shortly extend the discussion to sites with vector 
space of finite dimension d > 2, see also appendix [A] . 
A basis of the vector space V®^ of the lattice system is 
given by 



\SiS2 ■■ - Sn 



> = (5lr(4r---(4)^"ioo---o). (12) 



Recall that fermionic operators can be expressed in terms 
of Pauli matrices {(7^,(7^,(7^} by means of a Jordan- 
Wigner transformation. 



Cfc 



\k'<k 



at + i(Tl 



(13) 



The fermionic lattice system C is further characterized 
by a local fermionic Hamiltonian H . This is a Hamilto- 
nian that, when expressed in terms of the fermionic op- 
erators {cfc}, decomposes as the sum of terms involving 
only a small number of sites. As in Eq. ([l]), let j^f) S V**^ 
be a pure state of lattice C, 



m = 



E 



|SlS2 ■ 



■ Sn) 



(14) 



Here we will assume that \'^) is somehow related to the 
fermionic Hamiltonian H, for instance it is its ground 



state. Once more, we would like to find a varia- 
tional ansatz depending on 0{N) parameters to efh- 
ciently encode the tensor ^ containing the coefficients 
'i>siS2-siv of a pure state \^). 

One possibility would be to use a PEPS exactly as in 
the bosonic case. However, this might not be a good idea. 
Remember that the label fc € {1, 2, • • • , N} provides an 
order to the set of sites of C, whose position in the lattice 
is given by r = {x{k),y{k)). Two nearest neighbor sites 
fi and f2 on the square lattice might correspond to val- 
ues fci and ^2 that are far apart. Then, when expressed 
in terms of Pauli matrices, the local fermionic Hamilto- 
nian H will no longer look local. For instance, a nearest 
neighbor hopping term 



n 



(15) 



, fel<fe'<fe2 



develops a string of ir^'s. This might be harmful in two 
ways. On the one hand, the presence of strings of ci^'s 
would require important modifications in the algorithms 
to approximate the ground state of H with a PEPS, ei- 
ther through imaginary-time evolution or energy min- 
imization. On the other hand, it is unclear that the 
PEPS itself, which was originally designed as an ansatz 
for ground states of local bosonic Hamiltonians, will be 
as good an ansatz also for ground states of fermionic 
Hamiltonians, given that the latter are non-local when 
expressed in bosonic variables. 

Below we will explain how to modify the PEPS so 
that it is suitable to study fermionic systems (see also 
Refs. llOtT^ . Before, however, we introduce the notation 
necessary to deal with local vector spaces V of dimension 
d>2. 



E. Parity 

Fermionic systems are governed by Hamiltonians that 
preserve the parity of the fermionic particle number, to 
which we refer simply as ^parity'. That is, fermions can 
only be created or annihilated in pairs, and parity is a 
constant of motion. As a result, we can assume that 
the pure state I*) € V®^ of lattice £ has a well-defined 
parity, and observables 6 and reduced density matrices 
are block diagonal in parity. 

Let us consider again the vector space V of a single 
site, now with finite dimension d > 2. It is natural to 
decompose V as the direct sum of an even parity subspace 
V^+^ and an odd parity subspace V*^^\ 



V(+) 



(16) 



and to choose a basis of vectors with well-defined parity. 
Accordingly, the physical index s describing one such ba- 
sis is decomposed as s = {p,ap), where p € {—1,-1-1} is 
the parity and ap (denoted q;+ and a_) enumerates the 
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FIG. 5: (Color online) A tensor network is 'fermiomzed' by 
following two rules (see text): (a) all tensors in the network 
are chosen to be parity preserving (P is the parity operator 
acting on the indices of tensor T); (b) line crossings are re- 
placed with a special gate, a fermionic swap gate X . 



different basis states with parity p. The parity operator 
P then acts on this basis as 

P\p,ap) = p\p,ap). (17) 

In the case of spinless fermions, with a local dimension 
d = 2, the two possible states of a site are the local 
vacuum |0), signaling the absence of a fcrmion, and the 
state |1), corresponding to one fermion. In our notation 
these states read 

1(1,1)) ^|0), |(-1,1))^|1). (18) 

In the case of the t — J model there are three possible 
states per site, {|0), | f), | i)} S V, where | f) and | i) 
denote an electron with spin up and spin down, respec- 
tively. In our notation these states read 

1(1,1)) ^|0), |(-l,l))^|t>, |(-1,2))^|;). (19) 

Finally, in the case of the Hubbard model with local 
dimension d — 4 there is an additional state, | ti) = 
c|c| |0), corresponding to a doubly occupied site, which 
in our notation reads 

1(1,2))^ in). (20) 

In analogy with the physical index s — {p,ap), we 
introduce a parity operator also on the bond indices of 
a PEPS. Accordingly, the tensor A^^^^^ has bond indices 
u = {p, Qfp), etc. This means that e.g. bond index u can 
take values 

ue{(l,i),---(l,i?+),(-l,l),--- ,(-l,i?_)}, (21) 

where the bond dimension D is given hy D = + 
D_. The actual values of and Z3_ can be chosen at 
convenience. 



F. Fermionization rules 

Given a PEPS for bosonic systems, cf. Fig. [ija), in 
this work we obtain a PEPS for fermionic systems by 
applying the two rules used in Ref. [T2| to fermionize the 
MERA. These rules are applied both to the PEPS and 
to all related tensor networks that are involved e.g. in 
optimizing the ansatz or computing expectation values 
of local observables. 

Rule 1: Each tensor T in a tensor network is chosen 
to be parity preserving, i.e. 

Tiii2...JA/ = ifp(ji)p(j2)---p(H/) = -1, (22) 

where p{ik) G denotes the parity of the basis 

state labelled by z^, see Fig.[5ja). 

Rule 2: Each crossing of lines in the tensor network 
is replaced with a fermionic swap gate X, see Fig. [sj^b). 
This gate implements a fermionic exchange and has the 
form 

^«2iljlj2 = '^il,i2^i2 Jl'5'(*l7 *2), (23) 

with S'(ii,«2) given by 

Accordingly, starting from a PEPS for a bosonic sys- 
tem, a PEPS for a fermionic system is built as follows: 
(i) choose all the PEPS tensors {yll'^} to be parity in- 
variant. For instance, in the case of a bulk tensor A^^f^j^^, 
choose 

41s<i. = 0, if p{u}p{l)p{s)p{d}p{r) = -1; (25) 

and (a) introduce a fermionic swap gate X on any cross- 
ing of lines, as illustrated in Fig. [6] for a 3 x 3 lattice. 

Rule 1 is very convenient from a computational per- 
spective: it ensures that the parity of the wave function 
is exactly preserved during (otherwise approximate) cal- 
culations, while the block diagonal structure of tensors 
can be exploited to reduce computational costs. In addi- 
tion, Rule 1 is important in order to account for the an- 
tisymmetric character of fermionic wavefunctions (Rule 
2) in a simple way. However, we emphasize that tensor 
networks made of parity preserving tensors are also use- 
ful to describe bosonic systems (e.g., a Z2 invariant spin 
system, such as the quantum Ising model) and therefore 
Rule 1 is not what turns a bosonic tensor network ansatz 
into a fermionic one. 

Rule 2 accounts for the fermionic character of the ten- 
sor network, in the sense that it is employed to mimic the 
effect of anticommutators in a network of fermionic op- 
erators, as justified in Ref. ^'m the context of the MERA 
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(a) 





FIG. 6: (Color online) (a) Fermionic PEPS for the state 
I*) G V®" of a 3x 3 lattice £, obtained from the bosonic PEPS 
in Fig. [l][a) by replacing line crossings with fermionic swap 
gates. The coefHcients ^'^1^2 ■■'9 ^.re obtained by contract- 
ing this fermionic PEPS using standard tensor multiplication 
techniques (see text). The presence of fermionic swap gates 
introduces a complex structure of minus signs in ^'iiij - sg- 
(b) Bulk tensor A^^, which is chosen to be parity preserv- 
ing, Eq. (25 1, (c) Hermitian conjug ate of A^^. (d) Scalar 
product {^]^) written as a tensor network of reduced ten- 
sors a''^ defined in (e). Notice that, in a L x L lattice, O(L^) 
fermionic swap gates pervade the PEPS and yet, thanks to the 
jump move, Fig. |10| it is possible to write the scalar product 
('l/j^') in terms of only 0{L^) fermionic swap gates and in 
such a way that all of them are near a pair 

(^n^^nt) and 

can therefore be absorbed into the definition of the reduced 
tensor a^^ . [For a detailed derivation, replace line crossings 
with fermionic swap gates in Fig. [2]. As a result, there are 
no fermionic swap gates left in £, and therefore this tensor 
network can be contracted using the techniques employed for 
bosonic PEPS. 



(see also the note on terminology in the Introduction of 
this paper). 

Several additional remarks concerning fermionic ten- 
sor networks and their manipulations are in order. We 
start with a number of comments on fermionic tensor 
network representations that are relevant to the present 
formulation of the fermionic PEPS ansatz. 

(i) Fermionic order. The label k G {1,2,- •• , N} for 






FIG. 7: (Color onhne) (a) The fermionic PEPS in Fig.[6]corre- 
sponds to a specific choice of fermionic order, which appears in 



the definition of the local basis of Eq. ( 12 1 and Jordan- Wigner 
transformation of Eq. (13 1, (b) Another possible fermionic 
PEPS, associated to another fermionic order. The two tensor 
networks are not equivalent in that they cannot be mapped 
into each other by jump moves alone (where one is not al- 
lowed to drag a line over the end of another open line), but 
the mapping is possible if we allow some additional fermionic 
exchanges to modify the order of the open lines. Importantly, 
it can be seen that the tensors {A'*^} in both PEPS appear 
in the same way in any expectation value. In particular, they 
would be optimized using exactly the same figure of merits. 
This shows that the tensors {yl''^} are independent of the 
choice of fermionic order. 



the sites of C establishes an order on these sites. This 
order has been used to define a local basis of the Hilbert 



space in Eq. (12), and can also be used to translate the 



local fermionic lattice model into a non-local bosonic one 



through the Jordan- Wigner transformation of Eq. ( 13 ) 



(although this is not the strategy that we follow here) 
According to our prescription to graphically represent 
tensors and tensor networks, this order is also the order 
in which the open indices are drawn in Fig. |6][a). Notice 
that the structure of line crossings in the PEPS depends 
on the fermionic order and therefore the number and lo- 
cation of fermionic swap gates, see Fig. [7| Notice also 
that, in contrast with Refs. 1101131 we do not explicitly 
introduce fermionic operators (and corresponding order) 
on the bond indices of the PEPS, but use instead a simple 
graphical notation and two rules that already account for 
the complex pattern of fermionic-exchange minus signs in 
tensor 

(ii) Local operators. A local fermionic operator 6 is 
characterized by a tensor of components that describes 
the action of 6 on a given basis of states. For instance, 
in the simple case where each site has a vector space V 
of dimension d ~ 2, we expand a two-site operator 6 as 



Jjlj2)(«li2| 



(26) 
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where 



b-ij2>(*i*2| ^ (a^yU^^n0i02)(0i02|(S)'M«)^S (27) 

and d and b are the annihilation operators acting on the 
two sites (see appendix [A| for the d> 2 case). The coef- 



ficients Oi. 



are given by 



For instance, a hopping term d^b reads 

a^&= |li02>(0il2|, 
since the only non-vanishing coefficient is 

ooioi - (Ii02|a^6|0il2) 

= (0i02|aa^6 6^|0i02) = 1, 



(28) 
,(29) 

(30) 

(31) 
(32) 



which is obtained using the anticommutation relations. 
When the two-site operator o acts on two sites fci, /c2 G -C, 
the state j^*) of the system is modified into some other 
state I'i'') = o|\I'). This is implemented simply by con- 
necting the tensor (or tensor network) that represents 
^sjS2...s„ with the four-index tensor Og^^sk-^s'^ s'^ as in- 
dicated in Fig. [Sf^a). In particular, when the two sites 
are not contiguous in the fermionic order, |fc2 — > 1, 
a number of line crossings appear. This is reflected in 
the computation of the expectation value (^|6|4'), see 
Fig. |8];b)-(d). 

(iiij Parity changing tensors. Parity preserving tensors 
allow us to represent both states \'^) with even fermionic 
particle number (i.e. parity p = 1) and local operators 
6 that are parity preserving. But they also allow us to 
represent states with an odd fermionic particle number 
(parity p — —1) and parity changing operators. This is 
so because a parity changing tensor T, 



= if p(ii)p(i2) . . .p(«m) = 1, 



can be represented as a parity preserving tensor T, 



T, 



T 

J- 7 1 



(33) 



(34) 



where the additional index j only takes one value, j — 
{p,aP) = ( — 1,1). For instance, in order to represent a 
state 1^) e with an odd number of particles by 

means of a PEPS, an additional index j — (—1, 1) is 
attached to one of the PEPS tensors, see Fig.|9];a). Sim- 
ilarly, a fermionic annihilation operator c, which changes 
parity, can be represented as a tensor with three in- 
dices, one of which is fixed to j = (p, a^) = (—1,1), 
see Fig. ^h). It is sometimes computationally conve- 
nient (e.g. in the computation of two-point correlators. 
Fig. 16 1 to represent a parity preserving operator that 



is the product of two parity changing operators, such as 



the hopping operator d^b in Eq. (30), by two parity pre- 
serving tensors connected by a fixed index j = ( — 1,1), 
see Fig. ^h). 




FIG. 8: (Color online) (a) The state |*') = o|*) is described 
by a tensor of coefficients obtained by connecting tensors 
\[' and o. The example corresponds to a 3 x 4 lattice and 
a choice of fermionic order such that the two sites on which 
6 acts, which are nearest neighbors in jC, are not contigu- 
ous in the fermionic order (they occupy positions fci = 5 and 
^2 = 8). The line crossings (fermionic swap gates) that ap- 
pear in connecting tensors ^' and o can be interpreted as 
changes in the fermionic order needed in order to bring the 
two sites on which 6 acts together, then bringing them back 
to their original position, (b) Expectation value (4'|6|*I'). (c) 
The same expectation value in terms of a fermionic PEPS, 
(d) After some manipulation, an approximation of the expec- 
tation value is expressed in terms of an approximate two-site 
environment 5'''^''"' and tensors A^^^^, A^^^^\ A^^^\ A^^^^'' and 
o, as it was done for a bosonic system in Fig. [4] The only 
difference here is the presence of a few fermionic swap gates. 



(iv) Simplification. In the particular case of a pair 
of crossing lines i and j where index j only takes one 
possible value p of the parity, the fermionic swap gate X 
reduces to a product of two gates. Namely, to a product 
of two identity operators I® / if p = -1-1, and to a product 
of the parity P on line i and identity / on line j if p = — 1. 
It follows, for instance, that the jump move applied to a 
line i and a parity preserving tensor T with an index 
j — {p,a^) — (—1,1) (used e.g. to represent a parity 



changing tensor T in Eqs. (|33)-34| allows us to ignore 
index j in T at the price of applying the parity operator 
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FIG. 9: (Color online) (a) PEPS for a state |*), of a 3 x 3 
lattice, with odd parity p — —1. All PEPS tensors are par- 
ity preserving. One of the tensors has an additional index 
j = (—1,1) with fixed parity p = —1. (b) The tensor o of 
coefficients Oi^i^j^j^ for a two-site operator such as a hopping 
term a'h in Eq. ( |30[ ), that is the product of two parity chang- 
ing operators, and h, can be decomposed as the product of 
two parity-preserving tensors o' and o" connected by an index 
j = (—1, 1) with fixed parity p = —1. This decomposition is 
useful e.g. in the computation of a correlator involving two 
distant sites, see Fig. |16| Notice that parity preserving tensor 
o' can also be used independently of o" in order to represent 
the parity changing operator a. 




FIG. 10: (Color online) (a) As a result of rules 1 and 2, lines 
can jump over (parity invariant) tensors, with the fermionic 
swap gates traveling with the line crossings, (b) A parity 
changing tensor T is represented with a parity preserving ten- 
sor T with an extra index j = (—1,1) that remains open. In 
this case the jump rule introduces an additional fermionic 
swap gate X, involving j and the jumping line, that reduces 
to a parity operator P on the jumping line (see text), (c) A 
line can also be dragged over a set of open lines provided that 
the latter have even overall parity (in other words, provided 
they could be connected to a parity preserving tensor) 



P to line z, see Fig. 10 'b). This simplification appears to 



be useful e.g. in th e ca lculation of the expectation value 
(vE'|a'f6|*), see Fig. 



16 



A fermionic tensor network can be manipulated sim- 
ply by performing a sequence of tensor multiplications. 
It turns out, however, that by considering a special prop- 
erty (jump move) of fermionic tensor networks obeying 
Rules 1 and 2, the fermionic swap gates can be treated 
in a very special and advantageous way: they can be ig- 
nored until they correspond to a crossing of two indices 
connected to the same tensor, in which case they can 
be absorbed into the tensor using a low cost operation. 
As a result, in a fermionic tensor network algorithm we 
can follow the same sequence of tensor multiplications 
that we would have employed to manipulate its bosonic 
counterpart. In particular, we can use the same optimal 
sequence of contractions, with the same computational 
cost (to leading order). Let us discuss all these aspects 
in more detail. 

(i) Jump move. It follows from Rules 1 and 2 that in a 
fermionic tensor network, as in the bosonic case, lines can 



be dragged over tensors, see Fig. 10 a). This invariance 
under 'jump' moves allows us to modify fermionic tensor 
networks in such a way that the fermionic swap gates 
do not increase the leading cost of manipulations. For a 



proof of this property, which exploits the fact that ten- 
sors are parity-preserving, we refer again to Ref. [12] No- 
tice that the jump move does not include dragging a line 
over the open end of another line (this would amount to 
a change in the local fermionic order). The latter trans- 
formation is only allowed when it involves a set of open 
lines with even overall parity, as explained next. 

(ii) Open lines. By construction, line ends only appear 
in the diagrams grouped in such a way that they could 
be connected to a parity preserving tensor (that is, their 
combined parity is always even) without introducing ad- 
ditional crossings to the network. For instance, the set 
of all open legs of the PEPS in Figs, [ejja) or [gja) have 
even overall parity, since the PEPS as a whole has parity 
p — Another example is given by the open legs of an 
environment. Fig. |14| [c). Since such groups of open legs 
with even overall parity could be connected to a parity 
preserving tensor, the jump move extends naturally to 
them, as illustrated in Fig. 10 ^c). 



(iii) Absorption of a fermionic swap gate. Given a 



tensor T, with coefficients Ti-^...ij...ij^j , such that the two 
contiguous indices i,j are connected to a fermionic swap 
gate, it is possible to absorb the latter into the former 
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FIG. 11; (Color online) Multiplication of two tensors in a 
fermionic tensor network, (a) Example of two tensors Ti and 
T2 connected by two lines, (b) By means of a jump move, all 
lines standing between Ti and T2 can be dragged away, (c) 
Any fermionic swap gate involving the lines that connect Ti 
and T2 can be absorbed into e.g. tensor Ti, which is mapped 
into some Ti. (d) Finally tensors Ti and T2 are multiplied 
together, producing tensor T. The cost of all these manipu- 
lations is dominated by the last step. 



and produce a new tensor T with coefficients 



7^. ... — 71. 



(35) 



The cost of the absorption is thus just proportional to 
the number of coefficients in T and comparable to per- 
muting two adjacent indices in a bosonic tensor network. 
This cost is, in particular, smaller than the tensor mul- 
tiplication that may have produced T or a subsequent 
tensor multiplication involving T. Therefore the absorp- 
tion of fermionic swap gates only makes a sub-leading 
contribution to the total cost of contracting a fermionic 
tensor network (in the same way that the permutation 
of indices only makes a sub-leading contribution to the 
total cost of contracting a bosonic tensor network). 

(iv) Same leading cost. In a bosonic system, the con- 
traction of a tensor network is implemented by a sequence 
of multiplications of tensors, where each multiplication 
involves two tensors and the sequence is chosen carefully 
to minimize computational costs. In the fermionic case, 
it is natural to ask whether the presence of fermionic 
swap gates will modify the optimal sequence of tensor 
multiplications and increase the computational cost of 
the contraction. It turns out that one can always follow 
the same sequence of multiplications as in the bosonic 
case. Indeed, as illustrated in the example of Fig. [TT] in 
order to multiply two tensors, one can use jump moves 
and fermionic swap gate absorptions to eliminate any 
fermionic swap gate between the two tensors, which can 
then be multiplied together exactly as in the bosonic case. 
The leading cost of these manipulations is given by tensor 
multiplications, and it is therefore the same for bosonic 



and fermionic systems. 



G. Fermionic PEPS algorithms 

We have explained how to build a fermionic PEPS for 
the state j^") € of a fermionic lattice system. The 

next step is to consider algorithms to compute the ex- 
pectation value (6) of a local operator o from a fermionic 
PEPS, as well as to optimize the coefficients in its ten- 
sors {AI'^}. It turns out that these algorithms can be 
obtained by just introducing simple modifications to ex- 
isting bosonic PEPS algorithms. Here we discuss these 
modifications broadly. 

Let us consider first the computation of the scalar 
product (^l^*). As we did in the bosonic case, see Fig. [2] 
we build this scalar product by connecting the open in- 
dices of a PEPS for \-^) with those for a PEPS for (^-1 
obtained by Hermitian conjugation. The resulting ten- 
sor network may contain, in the fermionic CcLS6, cl large 
number of fermionic swap gates. However, an important 
point is that, by using the jump move of Fig. [lOf^a), the 
scalar product (^|^) can again be re-expressed in terms 
of a tensor network £ made of reduced tensors {a''^}, see 
Fig.gd). 

As shown in Fig. [6]Je), the fermionic reduced tensor 
differs from its bosonic counterpart by the presence of 
four fermionic swap gates. Importantly, in a finite or infi- 
nite lattice with open boundary conditions, all fermionic 
swap gates present in the scalar (4'|vE') are absorbed into 
the reduced tensor a'^^. Thus, the tensor network £ con- 
tains no fermionic swap gates. As a result, an approxi- 
mation to the scalar product (^'I'l') can be obtained by 
contracting £ with exactly the same approximate con- 
traction techniques (namely MPS techniques for finite 
PEPS^^ and either infinite MPS or CTM techniques for 
infinite PEP^i^ as in the bosonic case. [In some im- 
plementations, it is possible to lower computational costs 
by considering the components A^'^ and A'^^'^ of a'^^ sep- 
arately. In this case the approximate contraction tech- 
niques have to be slightly modified so as to account for 
the four fermionic swap gates involved in the definition 
of al*^]. 

Figure [12] considers a PEPS for a fermionic lattice sys- 
tem on a cylinder. Its scalar product can again be ex- 
pressed in terms of a tensor network £, made of reduced 
tensors {a''^}, that contains no line crossings. Therefore 
its contraction can also be performed exactly as in the 
bosonic case. Instead, in a fermionic lattice system on a 
torus, the analogous tensor network £ contains a num- 
ber of fermionic swap gates, see Fig. 13 However 



it 

can be seen that even in this case £ can still be coarse- 
grained using the same TERG techniques of Ref. 18 THl 
for the bosonic case, by properly coarse-graining the swap 
fermionic gates at eac h step of the coarse-graining. 

As discussed in Sec. II C given an environment £^^^'^^\ 
an approximation Q^^^'^^i can be used both for the com- 
putation of the expectation value (o) of a local operator 



12 




FIG. 12: (Color online) (a) PEPS for a state I*) of fermionic 
3x3 lattice with cylindrical boundary conditions, (b) Tensor 
network for the scalar product (c) Tensor network E 

of reduced tensors {a'*^}. As in a system with open boundary 
conditions, E has no line crossings and can be contracted as 
in the bosonic case. 




FIG. 13: (Color online) (a) PEPS for a state |*) of fermionic 
3x3 lattice with toroidal boundary conditions, (b) Tensor 
network for the scalar product (c) Tensor network E 

of reduced tensors {ffl''^}. On the torus, it is not possible to 
eliminate all line crossings in which must contain fermionic 
swap gates and theref ore di ffers from the bosonic case. How- 
ever, TERG techniquei^^'^ can still be applied by noting that 
the coarse-graining of the reduced tensors {a'*^} can take 
place independently of the fermionic swap gates, which are 
also coarse-grained in an obvious way. The (coarse-grained) 
fermionic swap gates only need to be absorbed into the rest 
of the tensor network at the latest coarse-graining step, when 
the torus has been reduced to only a few lattice sites. 



6 and for the optimization of the tensors {AI^} defining 
the PEPS. The approximate contraction of f [''i''^] lead- 
ing to is very similar to the contraction of E for 
the scalar product (^'|^'), and can again be accomplished 
in the fermionic case using the same techniques than in 
the bosonic case. 

From C/l''i''2lj obtaining an approximation to the ex- 
pectation value (\['|6|^') involves contracting the tensor 



network of Fig. [8[d) , which differs from its bosonic coun- 
terpart in the presence of 12 fermionic swap gates. An 
analogous tensor network is central to the update of the 
PEPS tensors {A^^} during an imaginary-time evolution 
towards the ground state |^gs) of a nearest neighbor 
fermionic Hamiltonian TJ, cf. Eq. ( [s]) . M ore details of 
this update will be provided in Sec. |IIl] in the context 
of an infinite lattice system, where also the computation 
of two-point correlators between distant sites will be ad- 
dressed. 

To summarize, the differences between bosonic and 
fermionic PEPS algorithms are in practice reduced to: (i) 
use of parity preserving tensors {AI'^} to efficiently en- 
code the state \^) G V®-^; and (ii) presence of fermionic 
swap gates in some of the tensor networks that need to be 
contracted, e.g. in order to compute the expected value 
of a local observable or optimize the tensors {^['^} of the 
variational ansatz. However, thanks to the jump move, 
the optimal sequence of tensor multiplications involved 
in the contraction of a given fermionic tensor network is 
the same than in the bosonic case, as is the incurred com- 
putational cost. All in all, we see that fermionic PEPS 
algorithms can be obtained by introducing a few mod- 
ifications to bosonic PEPS algorithms. This is further 
illustrated in the next section, where we provide more 



details on the specific PEPS algor ithm used in Sec. IV 
namely the iPEPS algorithm.fi32ni 



iPEPS ALGORITHM FOR INFINITE 
LATTICE SYSTEMS 



III. 



In Sec. |IV| the formulation of fermionic PEPS pre- 
sented in this paper is tested by computing the ground 
state of a number of models on an infinite lattice. In this 
section we provide additional details on the fermionic 
iPEPS ansatz and algorithm used for those compu- 
tations. We start by reviewing the bosonic iPEPS 
algorithm.liSsni Then we describe the modifications re- 
quired in order to address a fermionic system. 



A. iPEPS for bosonic systems 

The iPEPS ansatz exploits translation invariance of 
a system on an infinite lattice £ to store the state 
using only a small number of PEPS tensors, which are 
repeated throughout the lattice. Here we consider an 
infinite square lattice H and assume that the ground state 
of the system is invariant under translations by one site. 
We use an iPEPS made of copies of two tensors A and B 
that are distributed according to a checkerboard pattern, 
that is, 

j^[{x^x+2y)\ ^ j^^ ^[(x,:r+2a+l)] ^ X, y G Z. (36) 

The reason to use two different tensors A and i?, in- 
stead of a single tensor copied on all locations (as one 
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FIG. 14: (Color online) (a) Environment g^^^^^^ for two con- 
tiguous sites fi,f2 € jC, in terms of the reduced tensors a 
and b corresponding to the PEPS tensors A and B. (b) Ap- 
proximate environment C/I'^''^' expressed in terms of 10 ten- 
sors, corresponding t o the 10 shaded regions in (a), as used 
in t he C TM approactP^'^. Some of the tensor are used in 
Fig. 16 b). (c) Approximate environment g^^''-^^^ written in 
terms of 6 tensors {Ea}, as used in the computation of the 
expectation value of a two-site observable 6 (Fig. |16[ a)) and 
the update of tensors A and B during an imaginary-time evo- 
lution (Fig, [so]). 



would expect for a transl ation invariant ground state) is 
that the iPEPS algorithnli^^ requires that translation 
invariance be partially broken to the above checkerboard 
pattern during intermediate stages of the optimization of 
the ansatz. Invariance under translations by one site is 
(approximately) restored at the end of the optimization. 
Needless to say, the same ansatz is also valid for systems 
where the ground state has checkerboard order (invari- 
ance by diagonal shifts), as will be the case in some of 
the results in Sec. HVl 



B. Expectation values 

Let us first consider how to compute the expectation 
value of a two-site operator o, see Eq. ([2| . This requires 
computing both (*|*) and (*|o|\l>). 

The scalar product (^'l^') is expressed in terms of an 
infinite 2D tensor network £ of reduced tensors a and b 
distributed according to a checkerboard pattern. As in 
a finite system, an environment f [''i''^! for two nearest 
sites r*! , r2 G is then buil t from £ by removing two 



reduced tensors, see Fig. 14 |a)-(b). Contracting f*"! 



requires an approximation scheme, which produces an 
approximate environment Q^^^^^^ consisting of six tensors 
{Ea} connected through bond indices that take x differ- 
ent values, see Fig. [T4|^d ). Here x quantifies the degree 
of approximation in C/^^'"^!. One possible approximation 
scheme consists in using infinite MPS techniques, as was 
discussed in the original iPEPS algorithm in Ref. [T71 
Another possible approximation scheme consists in us- 
ing CTM techniques ,123 as discussed in the context of the 
iPEPS algorithm in Ref. 

An approximation (o)^ to the expectation value (6) of 
a two-site observable 6 is then obtained by computing 



(WO 



(37) 



where / denotes the identity operator on the space V (g) 
V of two sites of £ and we use the same approximate 
environment Q^^^''^^ to compute (vl/jolvl/)^ and {^>\i\'^')^ = 
(^1^)^, where the latter is computed as an expectation 

value for 6 — 1. The approximate value (o)^ will in 
general differ from the exact value (o). One expects, 
however, that in the limit of a large x o^is recovers the 
exact value. 



lim (o) 



(38) 



In practice, we compute 
uesof X, e.g. x & {10,20, ■ 



)^ for increasingly large val- 
, 100}, until the expectation 



value (6)^ does no longer depend substantially in x, and 
assume that this corresponds to (o). The cost of com- 
puting (6)^ using CTM techniques scales with the PEPS 
bond dimension D and the environment bond dimension 

X as OiD^x"")- 

Notice that for any value of the bond dimension D, a 
PEPS produces a variational energy (H) , 



(39) 



that is (H) > -Eoxact, where inexact is the exact ground 
state energy. Therefore, if we could compute (H) ex- 
actly, we would obtain an upper bound to -Ecxact, which 
could e.g. be compared with another upper bound cor- 
responding to another variational ansatz. However, the 
approximate value (H)^ is not guaranteed to be an upper 
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bound to i5„xacf In this paper we assume that {H)x is 
an upper bound to -Boxact once it has converged for large 
values of x- This assumption is seen to be correct in the 
case of free fermions in Sec. IIVI 



(b) 



a 





C. Simulation of imaginary-time evolution 



In the iPEPS algorithm of Refs. I17I20I an approxima- 
tion to the ground state is obtained by simulating an 
evolution in imaginary time, Eq. ([8|. The optimiza- 
tion of tensors A and B during the imaginary-time evo- 
lution is achieved by breaking the evolution into steps 
(Suzuki- Trotter decomposition) and then considering up- 
dates that involve a single link. There are four types 
of updates, corresponding to the four inequivalent links 
given by the bond indices m, I, d, r of tensor A. For each 
link, new tensors A' and B' are produced as a result of 
an optimization that aims to account for the action of 
a two-site gate g = exp{—h^^'^'^^^St) acting on that link. 
We refer to Refs. 1 171^ for a detailed explanation on how 
to use the approximate environment Cjl^i''^! to obtain up- 
dated tensors A' and B'. The algorithm proceeds by 
iteratively updating tensors A and B until they converge 
to some pair of tensors that depend on x and the time 
step St. A finite time step St introduces errors in the 
imaginary-time evolution. A better approximation to the 
ground state is obtained by gradually reducing St until 
its value does no longer affect significantly quantities of 
interest (e.g. the energy). For instance, this occurs for 
St — 10~^ in the simulations of Sec. llVl The cost of simu- 
lating an imaginary-time evolution is proportional to the 
number of time steps that are being simulated and scales 
as O(x'^-D^), since each update requires reconverging an 
approximate two-site environment Q^^'-'^^^ 

A simplified way of simulating an evolution in imagi- 
nary time was proposed in Ref. ,19, This simplified up- 
date does not involve the environment Q^^^^^^ and has a 
much lower cost per iteration, namely 0(D'^d^ + D'^cfi) 
when applied in a straightforward way and reducible to 
0{D'^(f + D^d!^ + D^d'^) after more careful considera- 
tions. On the one hand, ignoring the environment (that 
is, the rest of the wavefunction) implies that the update 
may not be optimal, and indeed there are cases, such 
as the 2D quantum Ising model near its critical point, 
where it produces a less accurate approximation to the 
ground state,'2Sl although we also found many cases where 
it only performs marginally worse than when using the 
environment. On the other hand, its much lower cost ac- 
celerates the simulations considerably. Notice, however, 
that the computation of expectation values with the re- 
sulting PEPS still requires computing Q^'^^'^^\ which has 
cost 0{x^D^). 



FIG. 15; (Color online) Reduced tensor a defined in terms of 
PEPS tensors A and A^ and four fermionic swap gates, (b) 
Reduced tensor 6 defined in terms of PEPS tensors B and 
and four fermionic swap gates. 





FIG. 16: (Color online) (a) Computation of an approximation 
{'I'|6|*I')^ to the expectation value (^'|6|^') of a local observ- 
able 6 acting on two contiguous sites. This tensor network 
only differs from the one in Fig. [ije) for a bosonic system in 
the presence of 12 fermionic swap gates, (b) Computation of 
two-point correlators. For concreteness, we consider the ex- 



pected value (cjcj) discussed in Eq. (41 1. The figure shows 



the tensor network representing (^'|c|cj|^)^, where sites i^j 
are in the same row of C but separated by 4 sites. Notice the 
line connecting the two tensors corresponding to operators c\ 
and Cj, which crosses a number of other lines introducing a 
number of fermionic swap gates. Since this line corresponds 
to an index j = (—1, 1) with well-defined parity p = —1, the 
fermionic swap gates simplify into I ® P. 



D. Fermionic iPEPS 

As in the bosonic case, the fermionic iPEPS ansatz ex- 
ploits translation invariance of a system on an infinite 
lattice C to store the state 1^*) using only a small num- 
ber of PEPS tensors, which are repeated throughout the 
lattice. Whether we are interested in approximating a 
ground state invariant under translations by one site or 
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with checkerboard order, we consider again just two (par- 
ity preserving) PEPS tensors A and B as in Eq. ([36]). 

It may not seem obvious that a fermionic iPEPS, char- 
acterized by tensors A and B, represents a state |^) of the 
infinite lattice C that is invariant under diagonal shifts 
(checkerboard order). Indeed, the presence of the ubiqui- 
tous fermionic swap gates, which are not homogeneously 
distributed across the tensor network corresponding to 
the iPEPS (see Fig. |6] for an illustration in the case of 
a finite PEPS) may seem incompatible with the checker- 
board order. However, the checkerboard order becomes 
manifest during the computation of expectation values 
for local observables. 

Given a fermionic iPEPS, the computation of the ex- 
pectation value of e.g. a two-site observable o is very 
similar to the bosonic case. Since there are no fermionic 
swap gates in the tensor network £ for the scalar prod- 
uct (^'l^') or the two site environment gl^^^^]^ these ten- 
sor networks can be contracted exactly in the same way 
as in a bosonic iPEPS. [If one wants to use the decom- 
position of the reduced tensors a and b in terms of A 



and B, Fig. 15 the presence of four fermionic swap gates 
needs to be taken into account]. In the simulations of 
Sec. |IV[ we have used the directional CTM approach 
discussed in Ref. [20! in order to produce an approximate 
environment C/I'"!'"^]^ Then an approximation (6)^ to the 
expected value (6) is computed using the tensor network 
of Fig. 16 'a). Figure 16 'b) describes the tensor network 



that needs to be contracted in order to compute a two- 
point correlators. 

The simulation of time evolution proceeds in a very 
similar way as in the bosonic iPEPS algorithm, with 
the difference that the tensor networks involved contain 
fermionic swap gates instead of simple line crossings. A 
detailed description of the two updates used in this paper 
is presented in appendix. |B] 

Before we move to presenting the results of ground 
state computations, we conclude this section with a sum- 
mary of the fermionic iPEPS algorithm: 

(i) Ansatz. The state of the system on an infi- 
nite square lattice C is encoded in two parity symmetric 
tensors A and B. These tensors depend on 0{dD'^) pa- 
rameters, where d is the dimension of the vector space 
V of one site of C and D is the bond dimensions of the 
iPEPS. 

(ii) Computation of expectation values. Given tensors 
A and B, the reduced tensors a and b are computed ac- 
cording to Fig. |15| From the reduced tensors a and 6, an 
approximation Q^^^^^^ to a two-site environment £^^^^^^ 



is obtained using the CTM algorithm, see Fig. 14 The 
computational cost scales as 0{x'^D^), where x is the 
bond dimension of the approximate environment. Then 
expectation values are computed by contracting small 
tensor networks involving Q^^^"^^^ . For instance, an ap- 
proximation (6)^ to the expected value (d), Eq. ( [37| , for 
an operator 6 on two nearest neighbor sites is obtained 
according to Fig. |16[a), whereas an approximation to a 



two-point correlator is obtained according to Fig. 16 'b). 

(iii) Approximation of the ground state. Starting from 
e.g. random tensors A and B, an imaginary-time evo- 
lution is used to find an approximation to the ground 
state of a local Hamiltonian H. The two updates of ap- 
pendix [B] can be used: the standard update, which in 
general produces better ground state approximations and 
has cost 0{x'^D^)', and the simplified update, which has 
a significantly lower cost and is often only marginally less 
accurate than the standard update. 



IV. RESULTS 

In this section we test the fermionic iPEPS algorithm, 
as summarized at the end of the last section, for sev- 
eral models of free and interacting fermions in an infinite 
square lattice £. We start by considering an exactly solv- 
able model of free spinless fermions, which allows us to 
compare the numerical results with the exact solution, 
and therefore assess the accuracy of the approach. It 
turns out that, similarly as in the 2D MERA,Q^ the ac- 
curacy of the numerical results depends on the amount 
of entanglement in the system. We also compare the 
use of the standard and simplified updates discussed in 
appendix [B] The second model, describing interacting 
spinless fermions on a square lattice, is no longer exactly 
solvable. We compare our results for the phase diagram 
with the Hartree-Fock (HF) solution from Ref. [29l and 
show that our results for large D = 6 are an improve- 
ment upon the HF result. The final example is the t ~ J 
model, where we compare energies obtained with iPEPS 
with previous variational Monte Carlo results based on 
Gutzwiller-projected ansatz wave functions. In all these 
examples we present iPEPS results for different bond di- 
mensions D, and study the convergence of the energy as 
a function of X; the bond dimension of the environment. 
In all these examples we chose = Z3_ = D/2, see Eq. 

pTl). 



A. Free spinless fermions 

1. Model 

The first model under consideration is an exactly solv- 
able model of free spinless fermions,!^ given by the 
Hamiltonian 



E 

(iJ) 



[4c,+H.c.-j{clc]+c,c,)]^2XY,4c^, (40) 



with (ij) denoting the sum over nearest-neighbor pairs, A 
the che mical potential, and 7 the pairing potential. Fig- 
ure 17 shows the exact phase diagram of the model. '^'^ For 



7 = the model reduces to the usual tight-binding model 
of free fermions, with a metal phase for < A < 2 and a 
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FIG. 17: (Color online) Phase diagram of the free-fermion 
model [40] as a function of chemical potential A and pairing 
potential 7. For 7 > 0, A > the system is superconducting. 
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FIG. 18: (Color online) Relative error of the ground state 
energy of the free-fermion model ( 40 1 as a function of A, for 



different values of 7 and D. The dimension x is 20 and 40 for 
D = 2 and D — i, respectively. 



band insulator for A > 2. Also the line A = 0, 7 > cor- 
responds to a metal phase with a one dimensional Fermi 
surface. For 7 > 0, A > the system is superconducting, 
with a critical phase for < A < 2 and a gapped phase 
for A > 2. 



2. iPEPS results 



Fig. 18 shows the relative error in the ground state en- 
ergy obtained by simulating an imaginary-time evolution 
with the standard update (cf. appendix [B|) . Results for 
bond dimensions D = 2 and D = A and at different lo- 
cations in the phase diagram are shown. In the gapped 
phase, A > 2, accurate results are already obtained for 
D = 2, and for Z? = 4 these accuracies are increased 
by one to three orders of magnitude. [A special case is 
the band insulator for 7 = 0, A > 2, (not shown) which 
corresponds to an unentangled or product ground state 
and can be reproduced exactly even with bond dimension 
D = 1.] The accuracy in the critical phase, A < 2, is of 
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FIG. 19: (Color online) Upper panels: Correlation function 
C(r-) = (clcj) as a function of distance r (in x-direction), see 
Eq. ( |41|, f or two different values of (7, A) of the free-fermion 
model ( |40| . Middle panels: Absolute value of the correlation 
function C(r) in semi-logarithmic scale. Lower panels: The 
difference between the simulation result C(r, D) and the exact 
result Cex{r) for different values of D. Notice that better ac- 
curacies are obtained in the gapped phase, which corresponds 
to a less entangled ground state. 



the order of 1% for D = 2. Increasing the bond dimension 
to D = 4, an order of magnitude is gained for 7 = 1, but 
the gain is smaller for 7 = 2. Finally, the free fermion 
regime with a ID Fermi surface (7 = 0, A < 2) is the 
most challenging case. Note that in this case, the entan- 
glement entropy exhibits a logarithmic correction to the 
area law.l^^MSI This shows that the accuracy of the results 
depends on the amount of entanglement in the sys tem, 
similarly to the findings with the 2D MERA.I21i2134] 

Next we study the accuracy obtained for the two-point 
correlator 



C{r) 



(41) 



where {x{j),y{j)) — {x{i) +r,y{i)), i.e. site j is in the 
same row as site j but separated by r — 1 columns, see 
Fig.[l6]^b). The iPEPS resuhs shown in Fig.[l9 



are seen 
. In the 



to approach the exact values with increasing D 
critical phase (left panels) only the correlations at short 
distances are reproduced accurately. The precision of the 
correlator at longer distance is rather poor for small D, 
but improves upon increasing D. In the gapped phase 
(right panels) the accuracy is clearly better, already for 
small bond dimension. For D = 6 accurate results are 
obtained up to distance r = 9 where the magnitude of 
the correlator is of the order 10^^. These results are ob- 
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FIG. 20: (Color online) Relative error of the ground state 
energy of the free-fermion model (401 as a function of x, the 
bond dimension of the environment, for D = 4. In these ex- 
amples the energy decreases monotonously with x, i-e. the en- 
ergy from the simulation is always larger than the true ground 
state energy. 
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FIG. 21: (Color online) Relative error of the ground state 
energy of the free-fermion model ( |40[ ) as a function of the 
number of renormalization steps in the CTM algorithm. The 
boundary tensors are initialized randomly, and the dimensions 
are _D = 4, and x = 32 (x = 16 for A = 1). The environment 
in the gapped phase converges considerably faster than in the 
critical phase. 



tained from simulations with the simplified update, and 
we checked that they are similar to the ones obtained 
with the standard update (for D = 2 and D — 4). We 
compare and discuss the two updates further below. 



3. Technical comments 

Since the iPEPS is a variational ansatz, the ground 
state energy (H) of a state represented by an iPEPS is 
an upper bound of the exact ground state energy E'exact- 
With increasing bond dimension D the true ground state 
can be represented more accurately, and consequently the 
energy (H) becomes lower, i.e. a better upper bound of 
E^^^^^. However, as discussed in the previous section, 
quantities such as the energy can only be extracted in 
an approximate way from the iPEPS. The error of this 
approximation depends on the bond dimension of the en- 
vironment X- In Fig- [20] the dependence of the relative 
error of the energy as a function of x is plotted for an 
iPEPS with bond dimension D = A. We observe that, 
with increasing the energy (i?)^ converges to some 
value that is indeed an upper bound for inexact- This is 
consistent with the assumption that {H)^^ has converged 
to the true energy {H) of our D ~ A iPEPS. In the other 
models analyzed in this section we will also assume that 
once {H)y. does no longer change with increasing Xi it 
has attained {H) and therefore is an upper bound to 
E^^aaf Fig. [20] shows that, in the present free model, the 
approximate energy (H)-^ monotonically decreases with 
increasing x, and therefore any value of {H)x is already 
an upper bound to inexact- However, this behavior is not 
true in general, as we will see further below. It is there- 
fore important to study the convergence in x in each case 
separately. 

Let us make a few comments about the convergence 
of the environment in the iPEPS algorithm, for fixed 
sets of tensors A and B which approximate the ground 



state of different phases. Figure 21 illustrates how the en- 
ergy converges with the number of renormalization steps 
in the directional CTM algorithrrP^ for different points 
in the phase diagram. The environment in the gapped 
phase clearly converges faster than in the critical phase. 
Here, the tensors at the boundary of the iPEPS have 
been chosen randomly at the beginning. Depending on 
the initial conditions of these boundary tensors the num- 
ber of steps needed to converge can vary significantly. 
In practice, when using the standard update, we use the 
environment from the previous imaginary time step as 
an initial condition, so that only a few renormalization 
steps are needed to re-converge the environment. Note 
also that for different initial boundary tensors the energy 
might converge to slightly different values (especially in 
the critical phase). It is therefore advisable to check re- 
sults for different initial boundary tensors. 

For this model, we have also compared the precision 
of the ground state energy obtained with the standard 
and simplified updates, see Fig. [22] The accuracies ob- 
tained with the standard update are typically slightly 
better than the ones obtained with simplified updates. 
However, the simulations with the simplified update are 
computationally considerably cheaper, because the en- 
vironment has to be computed only once at the end of 
the simulation for the evaluation of observables, whereas 
in the standard update the environment it has to be re- 
computed at each step of the imaginary-time evolution. 
The leading cost of the two methods is the same, but 
the computational cost differs by a rather large factor, 
which depends on the total number of imaginary-time 
steps. For the remaining examples of this paper we will 
consider only simplified updates. This allowed us to per- 
form simulations with D — 6 and x = 60 on a standard 
computer in roughly one day. We checked ioi D ~ 2 
and D — 4 that the results obtained with the two up- 
dates give similar accuracies. It is conceivable, however. 
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FIG. 22: (Color online) Comparison of the accuracies of the 
ground state energy of the free-fermion model ( 40 1 obtained 



with the two different updates described in appendix |B| 



that the difference in accuracy between the two update 
schemes become larger with increasing D, depending on 
the model under consideration. 




FIG. 23: (Color online) Phase diagram of the interacting spin- 
less fermion model ( |42[ ). V is the interaction strength, and 
n is the particle density (filling). The charge-density wave 
phase (checkerboard pattern) at exactly half-filling (n = 0.5) 
is separated from the metal phase by a first order phase transi- 
tion, with an intermediate region corresponding to phase sep- 
aration (PS). With increasing D, the boundary between the 
metal phase and the PS region moves away from the Hartree- 
Fock (HF) result from Ref. 1291 Dotted and dashed lines are 
a guide to the eye. The uncertainty of the phase boundaries 
obtained with iPEPS is smaller than the symbol size (in x- 
direction) . 



B. Interacting spinless fermions 

1. Model 

The second model under consideration is a model of in- 
teracting spinless fermions, which is not exactly solvable. 
It is defined by the Hamiltonian 



(ij) » (U> 

(42) 

with t ~ \ the nearest-neighbor hopping ampli- 
tude, V > Q the (repulsive) nearest-neighbor interaction 
strength, and /i the chemical potential. The Hartree-Fock 
(HF) phase diagram^^ as a function of V and particle 
density n is given Fig. [23] (solid line). The HF calculation 
predicts a gapped charge-density-wave (CDW) phase at 
half filling (n = 0.5) and a translational invariant normal 
state (metal phase) far away from half filling. In between 
these two phases they find a thermodynamically unsta- 
ble region, which we identify as a phase separation (PS) 
region, i.e. where the system splits into two parts, one in 
the metal phase and the other in the CDW phase. 



2. iPEPS results 

Our results for the phase diagram, obtained with the 
simplified update using an iPEPS with bond dimension 
D — A and I? = 6, are given by the squares and triangles 
in Fig. [23j respectively. The phase diagram obtained 
with iPEPS qualitatively agrees with the HF solution. 
However, with increasing D the phase boundary to PS 
region moves away from the HF result. In the following 



we explain how we computed the phase boundary and 
discuss the origin of this deviation. We focus only on the 
left half of the phase diagram, n < 0.5, since it is mirror 
symmetric with the n = 0.5 line. 

We determined the phase boundary between the metal 
phase and the PS region for several values of V, for = 4 
and D = 6. For each value of V, we studied the first 
order phase transition, which occurs at a certain value 
fi — II* . The iPEPS algorithm known to be particularly 
suitable to study first order phase tr ansiti ons, thanks to 
displaying some degree of hysteresis! ^^ * ^^ * Specifically, if 
we start a simulation with an iPEPS that represents a 
state in e.g. the metal phase, it will remain in the metal 
phase upon increasing ^, even for values (slightly) larger 
than n* (where the ground state is no longer metallic). 
This allows us to compute the energy of the metal phase 
in the region /i > /j,* even though the CDW-state has a 
lower energy, and vice versa. We can therefore compute 
the energies of the ground states of the two phases indi- 
vidually, and the phase transition occurs where the two 
energies cross, as shown in the upper panel of Fig. [24] for 
V = 2. At the transition point /i*, the density n jumps 
from a certain value n* < 0.5 in the metal phase to the 
value n — 0.5 in the CDW phase, as plotted in the lower 
panel of Fig. 



24 



For densities with values between n* 
and n = 0.5 there exists no homogenous ground state, 
i.e. the system exhibits phase separation (cf. Fig. 23 1. 



Let us now compare the energies obtained with iPEPS 



with the HF solution in Fig. 23 While the results in the 
CDW phase coincide (up to « 0.5%), we obtain lower 
energies for D = 6 than the HF result (of the order of 
3%). This leads to a shift of the transition point fi* to a 
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FIG. 25: (Color online) The convergence of the energy per site 
as a function of x is not monotonous. For large x the values 
do not seem to change significantly anymore. We associate 
an error bar to the values of the energy, depending on the 
convergence behavior in x- This error is smaller than the 
symbol sizes in Fig. |24[ 



FIG. 24: (Color online) Upper panel: Energy per site of the 
interacting fermion model ( 42 I as a function of chemical po- 
tentialu for V = 2, obtained by iPEPS and Hartree-Fock 
(HF)^ - The first order phase transition between the metal 
phase and the charge-density wave (CDW) phase occurs at a 
value fi* where the two corresponding energies cross. Lower 
panel: Particle density n as a function of chemical potential. 
At the first order phase transition point /i*, n jumps from a 
certain value n* in the metal phase to n = 0.5 in the CDW 
phase. For densities in between n* and n — 0.5 the system 
exhibits phase separation. The numbers in brackets indicate 
the uncertainty in the last digit. 



larger value of fj,, and thus the phase boundary appears 
at a larger value of n* than in the HF case. The error of 
the energy due to the finite x is smaller than the symbol 
size (of. Fig. 25 ). Within this error, we regard the results 
as "variational" energies, i.e. an upper bound of the true 
ground state energy. Accordingly, our D — 6 results in 
the metal phase are closer to the true ground state energy, 
and the phase boundary for D = 6 is an improvement 
upon the HF solution. 



3. Technical comments 

Note that the convergence of the energy with x iii 
Fig. [25] is not monotonous as in the previous examples 
in Fig. [20] Still, for large x the energies do not seem 
to change significantly anymore. For some simulations 
we observe that the environment does not converge to a 
fixed point, but rather oscillates slightly. In such a case 
the energy fluctuates around a certain value. We take 
this error also into account in our study. The final un- 
certainty of the phase boundary obtained by iPEPS is 
smaller than the symbol sizes (in x-direction) in Fig. 23 
All simulation results for this model were obtained with 



C. t — J model 

1. Model 

As a final example we consider the t — J model. 



+H.c.+j'^{SiSj--hihj)-fj, 



(43) 

with a = {tii} the spin index, fii — ^^c[„Cia the elec- 
tron density and Si the spin 1/2 operator on site i, and 



Cjcr — Cicj{l — (y^gCig). The t — J model is an effective 
model of the Hubbard model in the limit of strong on-site 
repulsion. The local Hilbert space of each site contains 
three basis states {|0), | f): I i-e. two electrons with 
opposite spins cannot occupy the same lattice site as in 
the Hubbard model. The t — J model is an important 
model in the context of high-Tc superconductivity and its 
phase diagram is still controversial. We focus here on the 
parameter tj J — Z which lies in the relevant parameter 
range of cuprate superconductors. At half filling, i.e. par- 
ticle density n — \^ the model corresponds to the antifer- 
romagnetic spin 1/2 Heisenberg model, where the ground 
state has long-range antiferromagnetic (Neel) order. Far 
away from half filling the model is a metal. At small, 
but finite doping x = 1 — n, several studies predict a 
dx^-y'i-^SMe superconducting phase, see e.g. Refs. 161371 - 
HUI and references therein. In studies with DMRG also 
a phase with stripe ordeil^^^^ was found at low doping. 
Some studies suggest that the formation of stripes is due 
to phase separation of the undoped antiferromagnet and 
the superconducting phase .'2SE3ES] 

Here we focus on the variational Monte Carlo (VMC) 
results from Ref. |3H1 based on Gutzwiller-projected 
ansatz wave functions. This study was done for finite 



20 



lattices of size 22 x 22, where the finite size corrections 
of the energy are estimated to be of the order of 10"'^ J. 
The Monte Carlo sampling error is of the same order 
of magnitude. At low doping (n > 0.9) the best varia- 
tional energies are obtained by an ansatz wave function 
including superconductivity and antiferromagnctic order. 
At larger doping 0.7<n<0.9a better variational en- 
ergy is obtained with an ansatz wave function that is su- 
perconducting without antiferromagnctic order. In the 
following we compare our energies with the best values 
obtained in this previous study. 
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Figure 26 shows the energy per site (with the chemical 
potential term subtracted) as a function of particle den- 
sity. One can see that the iPEPS results approach the 
VMC results from Ref. i38jwith increasing D. For D = 6 
the iPEPS energies are roughly 1% higher than the VMC 
energies, and for D = 8 the deviation is of the order of 
10~^ J, which is the same order of magnitude as the er- 
ror bar of the VMC results. Note also, that some of the 
energies for £) = 8 are lower than the VMC results. How- 
ever, the D = 8 results are not necessarily "variational" 
energies, as we discuss below. 

These preliminary results are encouraging for the fu- 
ture study of the phase diagram of the t — J model using 
PEPS algorithms, because the current energies with the 
largest dimension D = 8 are compatible with previous 
variational studies, and we expect to be able to increase 
D by properly exploiting the global symmetries of the 
model. We also would like to emphasize that iPEPS is 
a general ansatz, i.e. the same ansatz is used for any 
model on a square lattice, whereas in other variational 
studies the ansatz wave function is typically based on 
the specific physics of the model. It is remarkable that 
iPEPS, which starts from a random initial state, yields 
comparable energies as the ones obtained by specialized 
ansatz wave functions. Finally, it will be interesting to 
verify whether the phases obtained by iPEPS correspond 
to the ones predicted by the variational study, or whether 
a phase with different dominant correlations appears, 

e.g. str ipe -ordered phase, as previously found in DMRG 
studies.HSEll 



3. Technical comments 

Figure [27] shows that the energies as a function of x 
are sufficiently converged, both for = 4 and D ~ Q, 
i.e. the uncertainty due to a finite x is much smaller 
than the symbol sizes in Figj26] We can thus view the 
energies for I? < 6 as "variational" in the sense that 
they are an upper bound of the true ground state energy. 
However, for I? = 8 at finite doping the energy increases 
with increasing and it is at present not clear what 
value it will reach for larger (and presently unaffordable) 
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FIG. 26: (Color online) Upper panel; Energy per site in units 
of J of the t — J model as a function of filling n, with t/J — 3. 
With increasing D, the energies obtained by iPEPS are ap- 
proaching the values from the variational Monte Carlo (VMC) 
study in Ref. 1381 Note that the chemical potential term has 
been subtracted from the energy. The error bars of the VMC 
results are of the order 10~^ J. Lower panel: Relative devia- 
tion between the energies obtained by iPEPS and VMC. Full 
symbols indicate that iPEPS energy is lower than the VMC 
energy 



values X- Thus, in this case the energy is not believed to 
necessarily be an upper bound to the exact ground state 
energy. 



V. CONCLUSION 

In recent months several theoretical proposals have 
appeared describing fermionic versions of tensor net- 
work algorithms for 2D lattice s ystems, namely fermionic 
MERA^ and fermionic PEP^^SEMli, ^^^^ ^^^^^^ 

have explained how to obtain fermionic PEPS algorithms 
by applying the general fermionization procedure of ten- 
sor networks introduced in Ref. |T2] A highlight of our 
formulation of fermionic PEPS is its simplicity: it re- 
places the complexity involved in dealing with a network 
of fermionic operators with a tensor network built by fol- 
lowing two rules, namely the use of parity preserving ten- 
sors and the substitution of line crossings with fermionic 
swap gates. 
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FIG. 27: (Color online) Energy per site as a function of X- 
The energies for D < 6 do not seem to change significantly 
anymore for large Xi whereas for D — 8 the energy is still 
increasing at the largest value of x used. 



rate the internal global symmetries of a model (e.g. par- 
ticle number conservation) into the tensors and exploit 
their block structure to reduce computational costs.l^ 
This strategy is expected to be decisive for the charac- 
terization of the phase diagram of the t — J and Hubbard 
models. On the other hand, Monte Carlo sampling tech- 
niques could be used to reduce th e form al dependence 
of the simulation costs in D and ;^.li33Sl Finally, the use 
parallel computing on a large cluster would also lead to 
improved fermionic PEPS simulations. 
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Appendix A: Generalized fermionic operators 



We have then used the fermionic iPEPS algorithm to 
compute an approximation to the ground state of sev- 
eral fermionic models on an infinite lattice. By simu- 
lating an exactly solvable model of free fermions,'^ we 
have been able to see that even a fermionic PEPS with 
small bond dimension D {D = 2, 4) is already capable 
of reproducing up to several digits of the exact ground 
state energy, as well as two-point correlators. The results 
also showed that, similarl y to what had been observed 
with the fermionic MERAjl^IiS] the accuracy of the ap- 
proach depends on the amount of entanglement in the 
system. Generally speaking, gapped systems are less en- 
tangled than critical ones and, accordingly, a fermionic 
iPEPS with a given bond dimension D produces better 
accuracies for the former. The simulation of model of 
interacting spinless fermions has provided us with a first 
clear indication of the usefulness of the fermionic PEPS 
as a variational ansatz for the ground state of systems of 
interacting fermions. The fermionic PEPS approach re- 
produces the Hartree-Fock phase diagram^^, with metal 
and charge- density wave phases; but even with bond di- 
mension D = A and D = 6, it improves the ground state 
energies on the metallic phase and this results in a signif- 
icant shift of the phase boundary. Finally, results for the 
t — J model in the relevant parameter range for cuprate 
superconductors are particularly encouraging, given that 
fermionic PEPS, still at an early stage of development, 
already produce ground state energies comparable with 
those of previous variational studies .I^Sl 

The main limitation in present calculations is due to 
the scaling 0{x^D^) of simulation costs with the PEPS 
bond dimension D and the environment bond dimension 
X, which restricts D and x to relatively small values. 
There are, however, several ways in which larger values 
of D and x could become affordable, thereby leading to 
more accurate results. On the one hand, one can incorpo- 



In the case of a lattice where each site is described 
by a generic vector space V of finite dimension d, we 
decompose V into even and odd parity sectors, 

V9^v(+)®v(") (Al) 

and use an index s ~ {p, ap) to label a basis with well 
defined parity. 



P\p,ap) =p\p,ap), 



(A2) 



see Sec. 



II E We can also introduce a set {fs} of general- 
ized fermionic operators on each site, where fs is defined 
as 



fs^\0){s\, \s) = \{p,ap)). (A3) 

Then a local operator 6 acting on just one site can be 
expanded as 

6 = ^o,,,|^')(^l, k')(^l^/]'|0)(0|/.. (A4) 



Notice that fs is parity preserving if p{s) — +1 and 
parity changing if p(s) = — 1. Fermionic operators /ij^' 
and /is^' acting on two different sites rl, ^2 G £ fulfill 

=^(5i,S2)/i;^l/i?l, (A5) 

where 

-1 if p{si) = p{s2) = -1 



S{SI,S2) 



1 otherwise. 



(A6) 



A two-site operator 6 acting on sites fir2 G C, can then 
be written as 



0= XI Os2Sl4si|s'l4)(siS2| 



(A7) 



where 
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FIG. 28: (Color online) Definition of tensors A and B. 
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FIG. 29; (Color online) (a) Tensor A' is obtained from ten- 
sor A by crossing the legs s and d. (b) Tensor B' is obtained 
from tensor B by crossing the legs s and m. (c) Tensors Wa' 
and Ma« are obtained by joining the u, I, d indices and the 
r, s indices of A' , and doing a singular value decomposition 
of the resultant matrix. The singular values are included in 
the definition of tensor Ma* . (d) Tensors Nb and Wb are ob- 
tained by joining indices I, s and indices u, r,d oi B together, 
then performing a singular value decomposition of the result- 
ing matrix, then splitting the composed indices back apart. 
The singular values are included in tensor Nb- (e) Tensors 
W-^ and A/-4 are obtained by joining indices u, I, d and indices 
r,s of A, then performing a singular value decomposition of 
the resulting matrix, then splitting the composed indices back 
apart. The singular values are included in tensor iVfj. (f) 
Tensors Nb' and Ub' are obtained by joining indices s, I and 
indices u,r, d of B' , then performing a singular value decom- 
position of the resulting matrix, then splitting the composed 
indices back apart. The singular values are included in tensor 
Nb'. 




FIG. 30: (Color online) (a) Tensor Z is obtained by contract- 
ing a tensor network that contains the approximate environ- 
ment S''^""' = • • • , Se}, tensors Wa- , W^, Ub and Ub' , 
and a few fermionic swap gates, (b) Tensor Q is obtained from 
the contraction of Z, Ma' , Nb and the gate g. 



sen in order to best account for the action of the gate. 
Here we simply list the steps required in order to obtain 
the updated tensors A' and B' . For a justification of the 
schemes we refer to Refs. I14I17I201 We emphasize that 
the only difference between these updates and the ones 
used in a bosonic system is the presence of fermionic swap 
gates. In particular, if the fermionic swap gates are elimi- 
nated (by setting 5(ii,i2) = 1 in Eqs. (23)-24| we obtain 
update algorithms for bosonic PEPS. In some models. 



such as the quantum Ising model near criticality,'2Sl the 
standard update produces significantly more accurate re- 
sults than the simplified update, but this did not seem to 
be the case in the gapless phases studied in this paper. 
In all models away from criticality that we could test, the 
simplified update produces only marginally worse accura- 
cies. On the other hand, the much lower computational 
cost of the simplified update allows to consider larger 
bond dimension D than with the standard update. 

For concreteness, in the following we assume that the 
gate g is applied on a horizontal link where tensors A 
and B are at the left and right, respectively. Similar 
derivations apply to the other three types of links. 



Appendix B: Standard and simplified two-site 
horizontal update for fermionic gates 



1. Standard update 



This appendix describes how to update the tensors A 
and B that define a fermionic iPEPS. We consider the 
two strategies employed to obtain the results of Sec 



IV 



for 



namely (i) the standard update, used in Refs. I17|20 
bosonic systems, which requires the a ppr oximate envi- 



ronment = {El, - ■ ■ ,Ee} of Fig. [l4[d); and (ii) a 

simplified update, used in Ref.ll4[ which does not involve 
an environment. 

In these two schemes, a gate g — eyi'p{~h}^^^^^5t) is 
applied to the two sites ri,r2 G of a given link, with 
tensors A and B, and new tensors A' and B' are cho- 



Given tensors A and B, an approximate environment 
glrir-i] — {El, ■ ■ ■ ,Eq{ is obtained as described in Sec. 
|III[ first build the reduced tensors a and &, Fig. [15) which 
are the building blocks of the exact environment 
then use e.g. CTM technique to produce an approx- 
imation ^[^i^^l to f Pi*'^]^ 

On the other hand, it is convenient to introduce a num- 
ber of additional tensors. First compute tensors A and 
B according to Fig. 28 as well as tensors A' and B' ac- 
cording to Fig. 29'a,b). Then perform a singular value 
decomposition of tensors A' , B, A and B* as explained 
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vergence the process explained in the caption of Fig. 31 
Finally, the updated tensors A' and B' for the iPEPS 
are obtained as indicated in Fig. |32] 




FIG. 31: (Color online) The updated tensors M^. and N'g 
are obtained by iteratively solving two linearized equations. 
Specifically, we iterate the following two steps until conver- 
gence: (a) Given Nb and Nb' , we obtain M'^, by solving the 
linear equation M'j<^,R = S, where R and S are obtained as 
indicated in the diagram. Then, we set Ma' = M'^, , and 
compute = Mj^, . (b) Given Ma* and M^, we obtain 
N'b by solving the linear equation N'bR — S, where R and S 
are again obtained as indicated in the diagram. Then, we set 
Nb ~ N'b, and compute A'^^. = Ng. 




FIG. 32: (Color online) (a) Tensor A" is obtained from the 
contraction of Wa* and Af^. . (b) The updated tensor A' 
is computed by crossing the s and d indices of A" . (c) The 
updated tensor B' is obtained from the contraction of A'^^ and 
Ub. 



FIG. 33: (Color online) (a) Local detail of an iPEPS ex- 
pressed in terms of tensors Ta and Tb and four weight matri- 
ces Ai, . . . , A4, as required for the simplified update, (b)-(c) 
Relation with the usual iPEPS tensors A and B. 




\ 



FIG. 34: (Color online) (a) Tensor and the tensor network 
that defines it. (b) Tensors and Tb and weight matrix Ai, 
obtained from the singular value decomposition of B. 



in the caption of Fig. 29][c-f). Here we use the nota- 
tion T = WtDtUt = WtMt = NtUt for the singular 
value decomposition of a matrix T, where Wt and Ut 
are isometrics and Dt is the diagonal matrix of singu- 
lar values. Notice that the singular value decompositions 
are performed by regarding a tensor as a matrix after 
grouping its indices into two sets. 

Tensor Z is then obtained by contracting the tensor 
network of Fig. [30[ a) , which in turn is used to produce 
tensor Q in Fig.TsOfb). From tensors Z and Q, updated 
tensors M^. and N'^ are obtained by iterating until con- 



2. Simplified update 

For the simplified update from Ref. T^", the structure 
of the PEPS tensor network is slightly different to the 
one considered so far in this paper. Here the infinite 
PEPS is specified by two tensors F^i and Tb, and four 
diagonal matrices Ai, . . . , A4 with non-negative diagonal 
entries that assign weights to the indices of F^ and F b , 
see Fig. [33f^a). Notice that the usual expression in terms 
of tensors A and B can be recovered e.g. by multiplying 
the square root of the weight matrices Ai, . . . , A4 with the 
tensors F^ and Tb, see Fig. 33'b)-(c). 
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(b) 



r' 




FIG. 35: (Color online) (a) Tensor Tl is obtained by multi- 
plying with the inverses of A2, A3 and A4. (b) Tensor T'g 
is obtained by multiplying Tb with the inverses of A2, A3 and 
A4. 



The simplified update consists of the following steps: 
first, a tensor T\ is computed in analogous way as tensor 
A' in Fig. 29 |a). Then tensor 8 is obtained by contract- 



ing the network in Fig. 34 |a) , and subsequently decom- 



posed thro ugh a singular value decomposition as shown 
in Fig. 



34 



h). This results in tensors and F^ and the 
matrix of weights Ai, which corresponds to the singular 
values of Q. At this stage Ai is truncated into A'^, which 
keeps only the D largest diagonal entries of Ai. Then 
tensors F^ and Tb are also truncated accordingly. 

Next, tensors F^' and F^ are obtained from (the 
truncated version of) tensors F^ and F^ as shown in 
Fig. 35 a)- (b). Finally, tensor F^ is obtained from F^' 



again in an analogous way as A' in Fig. 32 b) 



The updated infinite PEPS is given in terms of the 
new tensors F^ and F^, and the set of weight matrices 
A'l, A2, A3 and A4. 
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